clip-02-05.jl

Load Julia packages (libraries) needed for the snippets in chapter 0

using StatisticalRethinking, Optim
#gr(size=(600,600));

snippet 3.2

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]

5-element Array{Float64,1}:
 0.43 
 0.43 
 0.722
 0.569
 0.865

snippet 3.3

Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);

10000-element Array{Float64,1}:
 0.35 
 0.708
 0.56 
 0.558
 0.816
 0.271
 0.616
 0.71 
 0.539
 0.448
 ⋮    
 0.589
 0.563
 0.829
 0.398
 0.422
 0.63 
 0.642
 0.467
 0.572

In StatisticalRethinkingJulia samples will always be stored in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

Object of type Chains, with data of type 10000×1×1 Array{Float64,3}

Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

2-element Array{ChainDataFrame,1}

Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean     │ std      │ naive_se   │ mcse       │ ess     │
│     │ Symbol     │ Float64  │ Float64  │ Float64    │ Float64    │ Any     │
├─────┼────────────┼──────────┼──────────┼────────────┼────────────┼─────────┤
│ 1   │ toss       │ 0.635981 │ 0.140213 │ 0.00140213 │ 0.00134517 │ 10000.0 │

Quantiles

│ Row │ parameters │ 2.5%    │ 25.0%   │ 50.0%   │ 75.0%   │ 97.5%   │
│     │ Symbol     │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼────────────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│ 1   │ toss       │ 0.341   │ 0.541   │ 0.643   │ 0.741   │ 0.878   │

Describe the chain

MCMCChains.describe(chn)

2-element Array{ChainDataFrame,1}

Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean     │ std      │ naive_se   │ mcse       │ ess     │
│     │ Symbol     │ Float64  │ Float64  │ Float64    │ Float64    │ Any     │
├─────┼────────────┼──────────┼──────────┼────────────┼────────────┼─────────┤
│ 1   │ toss       │ 0.635981 │ 0.140213 │ 0.00140213 │ 0.00134517 │ 10000.0 │

Quantiles

│ Row │ parameters │ 2.5%    │ 25.0%   │ 50.0%   │ 75.0%   │ 97.5%   │
│     │ Symbol     │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼────────────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│ 1   │ toss       │ 0.341   │ 0.541   │ 0.643   │ 0.741   │ 0.878   │

Plot the chain

plot(chn)

snippet 3.4

Create a vector to hold the plots so we can later combine them

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")